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NOMENCLATURE 


Greek Symbols: 

=1 for Fiber 

= 2 for Matrix 

= Small time interval 

= Small temperature difference 

fll / n 2 n 3 = Natural coordinates for the triangular 

finite element 

© = Angular coordinate 

0 = Field variable 

4 = Specific heat per unit volume 


Alphabetical Symbols: 


A 



D 

f 


[P] 

g 

h 

K 


B. 


K ,K ,K 
X Y z 


[k] 


N. 

1 


= Area of the 2-D finite element 
= Constants 
= Constants 

= Domain of the finite element 
= See table 2.1 
= Heat flux matrix 
=s See table 2.1 
= See table 2.1 
= Thermal conductivity 
= See table 2.1 
= Themal stiffness matrix 
s= Interpolating function 
= Radial coordinate 



= Radius of fiber 

r2 = Radius of equivalent cell 

3^,82 = Domain boundary segments 

T = Temperature 

Tq 5= Applied temperature 

t = Time 

t^ = Impulse duration 

X = Axial coordinate 


Superscripts: 

a , (a) = a - constituent 

(e) = Finite element (e) 


Subscripts : 


1 

2 

X 

r 

i / j / k 


=r Fiber 
= Matrix 

= In the axial direction 

= In the radial direction 

= Nodes of triangular finite element. 



ABSTRACT 


A formulation has been developed to study the 
phenomenon of one-dimensional Heat Conduction through 
Unidirectional/ Continuous Fiber Reinforced Composites. 
Further, a transient Finite Element procedure has been 
developed to solve the above formulation for variations 
in thermal properties and composite geometries. 

An initial boundary value problem has been 
selected and solved for temperature-dependent thenmal 
properties. System sensitivity to non-linearities aris- 
ing from temperature dependence of thermal properties 
has been examined extensively for composites having a 
constant Fiber»-Matrix volume fraction. 

The Finite Element Method results have been 
presented in the graphical form and discussed. It has 
been found and demonstrated that non-linear behaviour of 
thermal properties could significantly alter system 


responses 
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Chapter 1 


INTRODUCTION 


1.1 COMPOSITE MATERIALS 

A composite material, in general, could be defined, 
as one having two or more distinct constituent materials 
or phases. They offer outstanding mechanical properties 
coupled with flexibility in design and ease of fabrica- 
tion. The use of composite materials has become quite 
common, these days, in industries like aircraft, auto- 
mobile, sports-goods and electronics. Notable among the 
many composite materials, are the fiber composites which 
have the ability to make use of the outstanding strength, 
stiffness and low specific gravity of fibers such as glass, 
graphite or Kevlar. Also, in the case of fiber composites, 
the ability to synthesize the material structure during 
manufacture is a great advantage. 

1.1.1 Characteristics 


Composites consist of one or more discontinuous 
phases embedded in a continuous phase called the matrix. 
The discontinuous phase is usually harder, stronger and 
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primarily meant for reinforcement. An exception is the 
class of rubber-modified polymers, which have a rigid 
polymer matrix filled with rubber particles. 

The properties of composite materials are strongly 
influenced by the properties of their constituents, their 
distribution and the interaction between them. They also 
depend on the geometry of the reinforcement with reference 
to the system. Concentration, measured in terms of volume 
or weight fraction, is perhaps the parameter with greatest 
influence on composite properties. It being controllable 
during manufacture, desired properties of the composite 
could be created with great facility. The concentration 
distribution of the particles, which depends on their 
spatial relations, is a measure of the homogenity of the 
system. The orientation of the reinforcement affects the 
isotropy of the system, 

3.1.2 Classif ication 

The strengthening mechanism of a composite, 
strongly, depends on the geometry of the reinforcement 
and hence, it would be useful to classify them accordingly. 

Composite materials can be, in this light, 
broadly classified as fibrous composites and particulate 
composites. Fig. 1.1 gives a more detailed classifi- 
cation. 
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COMPOSITE MATERIALS 


J 

i - 

FIBER REINFORCED 
(FIBROUS COMPOSITES) 



SINGLE-LAYERED 
(Including composites 
having same orienta- 
tion and properties 
in each layer) 



CONTINUOUS-FIBER-REINPORCED 

» 


UNIDIRECTIONAL BIDIRECTIONAL 
REINFORCEMENT REINFORCEMENT 

(Woven reinfor- 
cement ) 


PARTICLE REINFORCED 
(PARTICULATE COMPOSITE) 

I 

RANDOM PREFERRED 

ORIENTATION ORIENTATION 

) 


LAMINATES HYBRIDS 


DISCONTINUOUS-F IBER-REINFORCED 



RANDOM PREFERRED 

ORIENTATION ORIENTATION 


MULTI-LAYERED 
(Angle Ply) 


F ig . 1.1 
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The chief characteristic of a particle is that 
it is non-fibrous in that, its length is not very much 
greater than its cross-sectional dimensions. As our study 
is related to fibrous composites, let us have a more 
detailed look at their characteristics, 

1 4 1 . 3 Fibrous Composites 

With normal engineering materials, imperfections 
in a direction normal to applied loads are particularly 
detrimental to their behaviour. Man-made fibres have 
considerably lesser flaws perpendicular to their lengths 
and are therefore, stronger, resulting in their use as 
reinforcements . 

The most important reinforcement fiber is E-glass 
because of its relatively lower cost. Boron, graphite 
and aramid polymer fibers (Kevlar 49) are useful for their 
high stiffness. 

Fibers, because of their small cross-sectional 
dimensions, are not directly utilized in engineering appli- 
cations but are used as fibrous composites by embedding 
them in matrix materials. The matrix serves to bind the 
fibers together, transfer loads to the fibers, and pro- 
tect them against environmental attack and damage due 


to handling 
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Fibrous composites can be broadly classified as 
single-layer and multi-layer (’angle-ply) composites. 

Single-laver " composites may/ actually, be made-up of 
several distinct layers, each layer having the same orien- 
tation and properties. The entire laminate, thus, behaves 
like a " Single-layer" composite. " Multi-layer" compo- 
sites, on the other hand, consists of several layers, 
each layer or lamina being a single-layer composite. The 
orientation of these single-layers is varied according 
to the desired design. 

When the constituent materials in each layer are 
the same, they are called laminates. Hybrid laminates 
are multi-layered composites with layers made up of diffe- 
rent materials. 

The chief advantages of employing fibrous compo- 
sites lie in their high strength-to-weight ratio and their 
controlled anisotropy which facilitates variation in 
desired ratios of properties in different directions. 

1.1,4 Fibrous Composites as Insulators 

The combination of high stiffness and strength 
with low weight enables fibrous composites to be natural 
candidates for aerospace applications. They are employed 
in supersonic transport to accommodate thermal stresses. 
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in sppce transport as heat shields and nose-tip insula- 
tors, and in space structures to reduce thermal deflec- 
tions . 

An outstanding example of their use as insulators 
is the case of the reusable surface insulation (rSI) tiles 
employed on the space shuttle orbiter [l]. Their cost 
effectiveness, resistance to thermal shoch and high insu- 
lative efficiency make them the best choice for such 
applications. The RSI tiles consist of silicon fibers 
dispersed in silicon carbide. 

Heat treated polyacrylonitrile (PAN) based carbon 
felt, in a rough laminar carbon matrix, gives a material 
with high resistance to thermal shock [ 2 ]. PAN-felt carbon, 
in chemically vapor deposited carbon matrices, also are 
good insulators [ 3 ]. 

1.1.5 Heat Diffusion in Fiber Reinforced Composites 

It has been observed that temperature and mois- 
ture are the primary environmental conditions that affect 
the behaviour of a structural composite during its service 
life. Hence, the analysis of thermal and/or moisture 
fields in a composite is of great importance. 

Conduction, parallel to the fibers of an uni- 
directional fibrous composite, is a problem of conside- 
rable practical importance for those cases in which such 
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composites are used as thermal insulators. Re-entry 
vehicle heat shields and nose-tips, in which the primary 
direction of heat flow is the fiber axis, are typical 
examples . 

The reduction of these problems to the one dimen- 
sional type entails greater computational efficiency. Most 
materials, used currently for these purposes, exhibit 
temperature dependent diffusive properties. The result- 
ing nonlinearities render analytical treatment intract- 
able. lA/hile considerable work has been carried out on 
the finite element analysis of conduction, the studies 
were limited to linear cases [4]. It is, therefore, of 
interest to analyse the non-linear heat transfer under 
such circumstances. 

finite element method in thermal analysis 

The Finite Element Method (PEM) was first applied 
to structural analysis. The approach is now recognised 
as applicable to any problem that can be formulated in 
the differential or variational form. This recognition 
is a boon to the analysis of thermal fields and coupled 
problems which arise in many applications. Aerospace 
structures, specially, can operate in severe thermal 
environs which significantly affect their structural 
design. Examples include convectively cooled structures 
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under development for hypersonic transport, thermal pro- 
tection systems for the space shuttle and large flexible 
space structures currently being designed for earth orbit. 
The design of these systems requires a strong interaction 
between thermal and structural analyses. The FEM offers 
the greatest potential for efficient coupling of the two 
[5]. 


1.2.1 Essence of the Method 

The heat diffusion equation could be written as, 


v^kv't + q_c|| = o (1.1) 

which is a classical example of a continuum mathematical 
problem. 


The numerical 'discretization' of such a conti- 
nuum problem with its boundary conditions 


T - T = 0 on ^ 


1 


Tc V T _ q = 0 on w. 


( 1 . 2 ) 


is the essence of all numerical solutions, where and 
portions of the boundary on which temperatures T 
or fluxes q are given and n represents the unit normal 
vector to the boundary. 


The finite element method can, in general, be 


defined as a process wherein. 
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(i) the continuous function T is approximated by 
a series of parameters and specified trial 
functions N^(x,y, 2 ) in the problem domain Q as 

T,wT sSN^a^j^j i = l...n (1.3) 

(ii) writing the differential equation (such as (1.1)) 
and its boundary conditions (such as (1.2)) in 
the general form: 

A(T) = 0 in Q 

(1.4) 

B (T) =0 in w 


The approximating equations, from which the 
solution is to be obtained, are written as a set 


W. A ( T ) dQ + S W. B ( T ) dto = 0 
£3 J to 1 


j = 1 n 


(1.5) 


where Wj and are suitable weighting functions which 
ensure that as n 


A ( T ) -* 
and g ( rp ) 

i.e., at all points, 
the exact solution. 


0 in Q 

( 1 . 6 ) 

0 in CO 

the approximating solution tends to 


The definite integral of equation (1.5) can be 
written as simply a sum of such integrals occurring on 
'subdomains' into which the whole 'domain' is divided 
(Fig. 1.2). Thus, if 




^ XEasMsas 
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m ^ Q 

Q = S Q and w = E w 

e=l e-1 

then, for all finite functions ( ) , we have, 

m 

/„ ( ) dQ = S S ( ) dQ 
y - G 

e=l Q 

m 

/ ()da)= E /g()da) 

e=l 0) 


( 1 .7 ) 


( 1 . 8 ) 


where, and u) ^ are associated with 'elements' into 

which we divide the problem. 


This allows the whole region to be divided into 
standard type of sub-regions, in which the parameters 
aj^ are usually the nodal values of the independent func- 
tion T and in which the tr-’'al functions are defined in a 
local manner. The integrals can, then, be evaluated 
element by element and the aoproximating eguations (such 
as (1.5)) obtained by a simple addition of element con- 


tributions . 
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1 . 3 PREVIEW OF PAST WORK 

On reviewing the work done in the area of con- 
duction, it was found that while literature is rich with 
various methods for the solution of linear problems, there 
is relatively little material available on non-linear 
problems . 

In 1951, Storm [s] showed that the mathematical 

condition for conversion of the 1-D non-linear conduction 

problem to the linear case is the constancy of [l/(KS)^'^^] 
dr "1. / 2 *1 

[log (S/K) J. The conversion, however, was found 
valid only for simple metals. 

Pattle [ 7 ], in 1959, presented a closed form 
solution to the non-linear problem using an integral 
method for the restrictive case of an internal heat source 
in an infinite medium whose thermal conductivity follows 
some power law K(t)'-^t’^. During the same year, Yang et al., 
[s], discussed the integral method for multi corrponent 
planar solids, 

Tsang [ 9 ], Andre-Talamon [lo] and Olsson [ll], 
in later years, discussed small perturbation methods to 
solve non-linear conduction problems. However, these 
methods are complicated to use, as it becomes difficult 
to obtain solutions to higher order perturbation terms. 
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Goodman [l2] and Imber [is] have described inte- 
gral or heat balance methods of considerably low accuracy. 
The latter amended the integral method to include phase 
change problems. 

Padovan [l4], in 1974 used complex series repre- 
sentations and properties of adjoint differential opera- 
tors to develop a finite element procedure for non-linear 
conduction, which is rather complicated to use. 

Leelamma et al., [ 4 ], in 1981, have pointed out 
the non-availability of finite element solutions to non- 
linear conduction problems, thus far, and proposed a per- 
turbation method for weakly non-linear cases of an isotro- 
pic material. 

Analysis of behaviour of composite materials 
began, rather, late. It was only in the early seventies 
(^Q 71 _X 973 ) that some research was carried out on diffu- 
sion in composites, Bedford and Stem [l5,l6] proposed 
a continuum theory for diffusion in composites, modelling 
the constituents as superimposed continue. The same 
theory was found effective for one-dimensional wave pro- 
pagation in composites in the direction of reinforcement. 
Hegemier and others [ 17,18] presented a continuum theory 
and an analvtical procedure for modelling wave propaga- 
tion in laminated and uni-directional fiber composites. 
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The theory takes care of mixture interaction and micro- 
structural information, but is restricted in application 
to wave propagation only. 

Horvay et al., [l9j/ first studied the case of 
transient heat conduction in laminated composites, 
omitting the case of fiber composites. Padovan [ 20 J 
gave an analytical solution for transient temperature 
distribution of a 2— D general anisotropic half space. 
Later, using a variational procedure, he developed a 
finite element approximation for steady heat conduotion 
in anisotropic media [21,22]. The methods, however, are 
not applicable for transient and non-linear situations. 

Haltman et al., [ 23 ], examined experimentally 
and numerically, steady heat conduction in a two phase, 
two-dimensional model of aluminium bars dispersed in 
foamed polyurethrane. H6=gemier et al., [24,25,26], gave 
various continuum theories related to fiber reinforced 
composites and laminated composites but, applicable to 
the phenomenon of wave propagation only. In 1975,Nayfeh 
[ 27 ] gave a continuum mixture theory for heat conduction 
in laminated wave guides, using approximations enabling 
analytical solutions. 

Chen et al . , [2s], calculated the thermal con- 
ductivity of composite materials containing thinly 
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dispersed spheroidal inclusions. The associated problems 
of heat flow across one such inclusion embedded in an 
infinite matrix were solved, in the presence of thermal 
boundary resistance. Maewal et al . , [29], developed a 
continuum model for wave propagation problems, and applied 
it for diffusion type processes also, in a laminated 
composite with periodic micro structure and obtained 
finite difference solutions. Again, Maewal et al., [30] 
and Murakami et al., [31] developed mixture theories for 
the process of heat transfer in unidirectional fibrous 
composites with periodic hexagonal microstructurek 

Patankar [32], gave a finite difference solution 
of heat conduction in the case of non-uniform thermal 
conductivity in a composite slab. Poon [33], reviews 
the various methods to solve problems of heat conduction 
in laminated composites. Murakami et al . , [ 34 -], extended 
their earlier work [30], to include non-linear material 
behaviour. Laura [ 35 ], gave analytical and numerical 
solutions to determine the temperature distribution in a 
composite rod with variable internal heat generation. 

It was in April 1983 , that Arun Aggarwal [ 36 ] 
developed a Finite Element Analysis package for the 
longitudinal heat diffusion in fiber-reinforced compo- 
sites. The work, however did not consider non-lineari- 
ties arising due to temperature dependent material 
properties . 
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The Present Work 

In the present work/ non-linearities arising due 
to temperature dependence of material properties have 
been taken care of as an extension of the work reported 
in [ 35 ]. An attempt has been made to show how non- 
linearities could affect the heat diffusion and tempera- 
ture distributions. For this purpose, the author deve- 
loped a Finite Element software package and executed it 
on the DEC 1090 system. 

An initial boundary value problem was selected 
and solved for non-linear material behaviour. Polynomials 
were used to represent the nature of non-linearity and 
the entire excercise was carried out for two different 
sets of materials. It was found that non-linearities 


could significantly alter the temperature distributions 



Chapter 2 


FORMULATION AND ANALYSIS 


2 . 1 THE FIELD PROBLEM 

A large class of field problems are governed by 
similar partial differential equations for the field 
variable 0. Heat conduction, wave propagation and con- 
vective diffusion are three examoles. 

Consider that it is desired to evaluate the 
field variable 0 in a three dimensional solution domain 
Q, bounded by surface S (Fig. ?.t). The field equation 
to be solved could, in general, be written as, 


A. (K 
dx X dx 



(K 

y ay 



(K 


M) 

z d z 


f (x,y,z) 

( 2 . 1 ) 


where K , K , K and f are given functions. The phy- 
!X! "y z 

sical interpretation of the parameters in the above equa' 
tion depends on the particular physical problem. Table 
(2.1) lists some typical problems and their parameters. 


The boundary conditions of (2.1) could be written 

as , 

0 = 0(x,y/z) on 


(2 . 2 ) 




Fig. 2>1 - 3-D SOlUTION DOMAW 



No. 


Problem 


F ield 
parameter 
0 


K 

> 

K 


K 


b 


1. Diffusion-flow 
in porous 
media 

2. Electric 
conduction 


Hydraulic 

head 


Voltage 


Hydraulic Internal Boundary 
conductivity source flow 
flow 


Electric 

conductivity 


Internal 

current 

source 


Exter- 

nally 

applied 

boundary 

current 


3. Electro-static 
field 


4. Heat 

conduction 


Electric 

force 

field 

intensity 

Tempera- 

ture 


Permitivity 


Internal 

current 

source 


Thermal Internal Boundary 

conductivity heat heat 

genera- genera- 
tion tion 


Convec- 

tive 

heat 

trans- 

fer 

coeffi- 

cient 


Torsion 


Seepage 


Stress 

function 


Pressure 


Reciprocal 
of shear 
modulus 


Permea- 

bility 


Angle of 
twist 
per unit 
length 

Internal 

flow 

source 


7 . Magnetostatics Magno- Magnetic 

motive permea- 

force bility 


Internal Exter- 
magnetic nally 


field 

source 


applied 

magne- 

tic 

field 

inten- 

sity 


TABLE 2.1 : I DENTIFICATION OF PHYSICAL PARAMETERS 
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and 


K 


X 




n 


X 


+ K, 


^ 0 


■y bY \ 


+ K 


2 d2 


n 

z 


+ g (x,Y,z) + h (x,v,z) = 0 on S 2 (2.3) 

where g and h are known apriori and 

direction cosines of the outward normal to the surface. 


2 .2 FORMULATION 
2.2.1 Composition 

Consider a two phase material whose phases are 
of cylindrical shape, with all generators oriented in 
the same direction. For convenience, let us concern 
ourselves with materials of cylindrical shape with gene- 
rators parallel to the phase region generators. It is 
also assumed that phase cylindrical regions are unin- 
terrupted along their lengths. 

In the general fiber reinforced material, there 
is no specific geometrical distinction between the two 
phases. If we consider one of the phases as consisting 
of cylinders embedded in another phase, we will call the 
embedded cylinders the FIBRES and the phase in which they 
are embedded, the MATRIX. 
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2.2.2 GeoTTie'bry 

If we assume the fibers to be of a specific shape 
(in this case circular, say), the geometry is one of regu- 
lar arrays of identical fibers in a matrix. The most 
general ^cinds of arrays are perhaps the square, rectan- 
gular and hexagonal arrays. 

For the present study, we consider a hexagonal 
array of circular, cylindrical fibers ( constituent- 1) in 
a matrix (constituent-2), as shown in Fig. 2.2. 

Using cylindrical polar coordinates, the domain 
could be described as occupying a space 0 <_ x ^ L , 
0^r<<» , O<0<2 r (piq. 2.3). 

2.2.3 Formulation 

We consider the temperature or heat flux boundary 
conditions at x = 0 and at x = L to be such that the tem- 
perature distribution is similar in each hexagonal cell. 

As a consequence, the external boundary of the unit cell 
becomes a line of symmetry, and the component of the heat 
flux normal to the boundary vanishes. It is further 
assumed that there is no thermal resistance across the 
fiber-matrix Interface. 

To describe the heat conduction in a typical 
unit cell, the hexagonal boundary is approximated by a 






^ MATRIX SHEL^ 


Fig. 72 - HEXAGONAL ARRAY 
COMPOSITE CYLWCER 





CONSTRUCTION 










'v 



1' 



circle 


2 4" 


. As a result, the temperature distribution within, 
the cell is axisymmetric with zero heat flux normal to 
the outer circular boundary of the matrix. 


The initial boundary value problem that des- 
cribes transient heat conduction in the cell is then 
given by the following equations: 


(a) Conse-rvation of energy 


1 

r 


A. 

dr 






dx 




(a) 


4 


(a) ^ 

dt 


(a) 


(2.4) 


(b) 


Constitutive equations 


r . K 

■X ~ ~ X <ix 

„ K 

r dr 


(2.5a) 
(2 .5b) 


(c) Interface 

conditions 


yd) 

(x,r-^ , t) = T^^^ (x,r^ , t) 

(2.6a) 


(x, r^ , t) = ^ 

( 2 . 6b ) 

( d ) Symmetry 

condition 



(x, r 2 / t) =0 

(2.7) 


(e) Initial and boundary data. 

. . (OC) (Of') ,,(0f.) m(a) 

The quantities q^ . g^ .4 

K and K in the above equations are, respectively, 

the heat flux components, heat capacity, temperature and 
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thermal conductivity components of the a ( = 1 or 2) 
constituent. 


The problem has been solved to determine the 
evolution of the temperature field in the half space x > 0 
subject to a step function at the boundary. 


T (O, r, t) = T^ , 0 £ t < t^ 

0 . t > tQ 


(2.3) 


2 . 3 FINITE ELEMENT FORMULJ^TION 

The governing equations for the axi symmetric 
heat conduction problem could be written as 


1 A (K r (K ^) = IL (2.9) 

? ^ r 3r^ + px ^ X dx ^ ^ dt 

along with its boundary conditions (Fig. 2.4). 

As we consider the temperature evolution in the 
half space x > 0 , the boundary conditions at the two ends 

become 


At X = 0 


T = T 


0 


0 <_ t < t^ 


( 2 . 10 ) 


0 


t > t, 


'O 


At X 


id 

5x 


0 


By symmetry of the composite cylinder assembl- 
age, the heat flux normal to the cylindrical surface of 
the unit cell vanishes. Thus, 
at r = r 2 ^ - 0 


( 2 . 11 ) 
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Due to the symtnetry of the geometry and boundary 
conditions, it would suffice, if analysis were carried out 
for one half of the cylindrical cross-section. Thus, due 
to symmetry , 

at r = 0 , ~ = 0 and T is finite (2.12) 

dr 

We shall resort to the Galerhin's approach to 
derive the FEM equations. 


2.3,1 Finite Elemen t Discretization 

without deciding about the type of element to 
use, we can write 

T*®* (x, r, t) = (x,rjf |T(t)J 

, It] 

where, the interpolating functions ISI^^ need to preserve 
C° continuity only i.e., continuity of the field variable 
at the element interfaces. 


The Galerkin's criterion could then be applied 


to obtain 


, r 1 d /T, ^ 

^(e) ^i ^ r eir r 


(e) 


r 


■) + (K 


'b T 


(e) 


3x X 6x 


^T 

St 


4 ^ ] dD^^' = O 


.(e) 


(2.14) 


,(e) . 


where i = l,2,...n and D^®' is the domain of element e . 


For axisymmetry/ 

dD^ =2 T r dA 


(2.15) 
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Using this in (2.14), we get, 

/; N. [ i ^ (K V- ^ 


D 


( e) 


i " r dr r" dr 


■) + ^ (K. 


(e) 


4 


St 


( e ) , 


St 


X dx 

( e) 


-] 2 " r = O 


(2 . 16) 


Integrating, by parts, the first term w.r.t. 'r* and the 
second term w.r.t. 'x', combining and rearranging the 
terms, one obtains 


,(e) 


,(e) 


;; [n. ^ (K^r ^ ) + N, ^ (K„r|^ )] dA^®^ 

( e) 


1 /r r dr 


D 


i,x X dx 
(e) 


+ JS 4 r N. li dA^®^ 


D 


( e) 


i dt 


X 


= J 


top 


N. (K r 


dT 


(e) 


bottom 


i r dr 


'right 


•left 


dx + 


^right 

+ J N, (K .r ) 

^left 


X X dx 


X 


top 


! X 


dr 


not tom 


or/ 


(2.17) 




D 


+ ^ ’^i lit 

d'®’ 




,(e) 


jzT (K_r 

.(e) 


I I ■ K 3r 1 ' ■ - x\ / u S 

i " r hr r X dx x" 


(2 .18) 
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Introducing the boundary condition for flux on the element 

,( e) 


boundaries ^2 ' ws havr 


K r 


^i-pCe) .^(e) 

„ K r ^ 

r X gx 


r 

Substituting in (2. IS), we get 

,(e) 


+ r q (t) = 0 for t > O 


(2 . 19 ) 


.(e). 


SS [ N. K r + N, ^ K,r dA^^^ + 


D 


J L ^ n j-— 


+ // H r 

D(e) 


i,x x-^ ax 

dA^®^ + jZT’ N. r q dS- = 0 
<,(e) ^ ^ 

^2 

(2.20) 


Using from (2.13) in (2.20), 

IS [k r N. + K r ih ;N,x|'^] dA^®^ + 

( ) r i,r ^ X , , , 

D ^ 

( 1 $T i 


+ IS [|1 r N. dA^®n [t] 

( e) 


D 


-- ( 0 ) 

+ J2r N. r q dS 2 = O 


( 2 . 21 ) 


2 

where T = ‘^ / the time derivative of the temperature, 
t: 

Casting (2.21) into the matrix form, we obtain, 

;/ (K^r (N,r| |N,r-'^ ■*■ 1 n,x| + 


D 


(e) 


+ IS 


D 


(e) 


(U r 


nU 


+ j/ (n - r q dS2 


(e) 


= 0 


( 2 . 22 ) 
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[K](e) ,p-;(e) 

This is the basic element matrix equation. 

For each element, 

the thermal stiffness matrix, 

= /; (K r ^N,r ! /N,ri^ + 

-.(e) ^ J i 


D 


+ K r }n,xI [n,xI^) dA 
^ < J I I 

the thermal capacitance matrix 


(e) 


u r iHi 'hF dA<«> 
(e) ' J ( ^ 

D^e; 


and the heat flux matrix 


^ jN] r q*®' d^®' 


e) 


= 0 (2.2 3) 


(2.24) 


(2.25) 


(2.26) 


Our problem is an axisymmetric case, which 
requires continuity of the field variable alone. Let us 
choose triangular ring elements (Figs. 2.5, 2.6) with 
linear variation of field variable T along the element 
boundaries as it satisfies both C° compatibility and 
completeness requirements. For the three node element 
we have. 


T 


(e) 


T 

('''I'i ! 

' T^(t) ' 

n2 V 

T2(t) 1 

[ n3 i * 

^T3(t)J 


(x,r, t) 


(2.27) 
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For linear variation of T along the element boundaries ^ 


>^1 




(r - r^. x.+(r^.x^ - r^x^. )+(r^x - r x^) 
)+(r^x^ - r^Xj)+(r^Xi- 


The denominator of is equal to 2A where A is 
the area of the triangle. Likewise, it is also the deno- 
minator of ri 2 ^ 3 * 

Hence, 


Ho = N. 


(r^x - r x^) + (r x^ - r^x) + (r^Xi - r^^x^) 


2A 


and 


Hj = N3 = 


ir^x^ - r.Xj^)+(r. X - r x^)H-(r x. - r. x ) 

(2.28) 


2A 


Also , 


We 


N 


N 


N 


N 


ll +^2 + ^3 

. Nj + Nj 

+ = 

1 

1 a consequence of the 

above 

def initions , 

,NNi _ 

2A 



1/r ar 

“ 2A 


f)^2 

(x^ - ’'1 > 

“ 2A 


^2 ,r “ r 

2A 



(Xj - Xj) 

®3 

(2 

^3,r dr 

2A 

" 2A 


, and 

are constants . 



( b - r • 
k J 

) 


^l,x doT" 

2A 

~ 2A 


(ri - 

— * = ^TT* 




2 ,x 


2A 


L ikewise 
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N 


3,x 


!±i-l 

2 a 


_3 
2 A 


(2.30) 


where and are constants - 

We now proceed to evaluate the elemental 

matrices . 


^ ^ 'thermal conductivity (stiffness) matrix : 

We have 

[k]***^^ = // (K r ^N,r? / N,rl'^ 

j3(e) i ) L J 

+ K r jN,xi (N,x\'^)dA^®^ (2.24) 

^ L J L j 

Our thermal conductivity K is temperature depen- 
dent for its value. As suggested by Abel and Desai [46], 
assuming piecewise linearity for K i.e., a constant but 
different thermal conductivity for each increment, we 
could proceed as though the case were linear. Using this 
ideology and substituting for N, r and N,x from (2.29) 
and ( 2 , 30 ) we obta in , 



(2.31) 
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Alq'6'tiraic manipulation will prove r = r. r)^ 4 -r.Y ]2 i* 
Hence, 


[k]<' 


K 


// , 


spc^ BjBj+C^Cj 8^63+ CjC3 


Symm, 


2 2 
b;+c 

2 


2 2 i 

®3'^^3 j 


'^1 + '^j ''2 + “"k ’’3’ 


( e) 


... c. (2.32) 

For a triangular element, the integration formula is 


, T j.(p) '’■■') 

S n -1 n n T) « QA — 

^(e) ^ ^ (a+3+r+2).’ 


2A 


(2.33) 


Employing this, (2.32) transforms to 

>( 


[k] 


(e) K 


(e) 


4A 


(e)' 


B?+c2 B,B3+CjC3' 


Symm. 


® 2'''^2 ® 2 ® 3 ‘*'^ 2^3 


2 2 

®3-^S J 

(e) 




D 


[l<] 


(e) K 


(e) 


(e) 


12A 


(e)' 


2 2 
®2+S 


(ri+rj+rj^) 


B 2 B 3 +C 2 CC 


t.2 „2 

®3’^^3 


(2.34) 
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In equations (2.26) and (2-28), has been 

used ns IF it were a constant and eaual in both axial and 
radial directions. This has been done for convenience of 
expression. During solution, it is necessary to club the 
value of radial conductivity with the 'S'* terms and axial 
conductivity with the ' C' terms. 


( i i ) Thermal capacitance matrix 
We have 

(e) 


= /; 4 r {n U d-A^®^ 

_(e) Ml] 


(2.25) 


D 


Extending the argument of piecewise linearity to the 
thermal capacitance and considering it constant over an 

element, we have 


[c] 


(e) ^ ^(e) 


jj ‘n h ^ h ”2 

“ . r. nOaA<-> 


// 


+ >21, n3 


-n, ^13 

n? 


(e) 


’’i + q '’2 ''3’'^* 

Employing the Integration formula, this reaucea to 
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^ 60 

( 6r^+2rj+2r^) (2r^+2rj+r^) (2r^+rj+2r^) 

( 2r. +6r .+2r, ) (r. +2r .+2r, ) 

1 J K 1 J K 

Symm. (2r^+2rj+6r^) 

..... (2 . 35 ) 


( i ii ) Element heat flux vector ; 
From (2.26), 



'2 


This term exists only for elements at the boun- 
clary where the heat flux is specified. Cbnsider that on 
the boundary edge the heat flux is constant for any ele 
ment. This assumption simplifies the analysis. 

Thus , 

4 


= (t) 

= q'®' (t) 




hi *^1 ’^+0 j 



( 2 . 36 ) 

For a linear element, the appropriate integration fomula 


is 
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a 



13' 

(fTf+P+1 )/ 


whoro is the length of the linear element. 

Using this, (?.36) transforms to. 


{f| 


= q<'^'(t) S. 


r 

-M. 


r^/3 + rj/6 

r,/6 + r ./3 


(2. 37) 


• - -r j.y 

where ^i--! i® length of the boundary edge of the 


element. 


2.3.2 General Matrix Equations 

The element matrices as given above are calculated 
for each element and assembled together to give the system 
of simultaneous linear equations in the matrix form as 

[k] [t| + [c] { T } + (f} = 0 (2.38) 

Our problem being a transient case, we are con~ 
fronted with time derivatives of temperatures. Zienkiwicz 
[ 37 ] shows that either finite difference aonroximations 
or finite elements in time (use a shape function defined 
in the time domain) may be used for computer solution 
and that the methods are equivalent. 

We will adopt the former method, replacing time 
derivatives by finite difference approximations. 

From ( 3 . 38 ) at time t, 
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[Kt] [t^I + [C^] 
At time t“5t 


( 2 , 39 ) 


[’'t-eJ K-atj + = ° 


( 2 . 40 ) 


Then, 

Ki = 

r 

Tt- 

] + “ 
6tj ^ 2 


Or, 

= 

2 

( $' T - 

^ I t) 

Pt-6t 


rt] = h ^ - i\~ 6 t] lK- 5 t] ^ 2 . 41 ) 

iSubsti huting eguation (2,41) in (2,39) we have, 

[K,] ^ [C,] (|^( 

H. = 0 

Or, 

+ 5:^ C^tj "" ^^t^ ^ ( '^t-6t} ?t |'^t-6t j ^ 

-[Pt\ (2,42) 

F rom ( 2 . 40 ) , 


t-6tJ t.\-6tj “■ |^t-6t} ^ 

(2.43) 


Substituting this in (2.42), we obtain, 

([k^] +1^ [cj) - [c^](-[c^.,t]''[K^-5t] it-stj 
- [=t-6th' (h-6t(> + h l^t-«t} - (,*'4 


Rearranging, 
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([K,] + [C^]) = ( |-[cj. 

- 

-w* 

For most materials the variation of thermal conductivity 
with temperature is much more pronounced than that of 
the thermal capacity. Considering, therefore, that this 
variation be neglected in comnarison, we can write 

= [c^l = [c] ( 2 .«) 

.'.([k^] + 1^ [c]) ■= [c] - 

- ( ptj Ih-St} > 

This equation is valid for transient analysis. 

For steady state problems, the time derivative of tempe- 
rature! is zero, reducing (2.38) to the form of the Fourier 
law of conduction, 

[k] T j + |f^ = 0 (2.47) 

Solution o(- this equation gives the steady state results 
in one step of computation. 

2.4 SOLUTION PACKAGE 

A computer program has been formulated to solve 
the matrix eauatlons In TORTEWN. The program enables 
holding of toundary nodes at constant temperatures. 
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As we arc solving for the temperature field/ it would 
be of interest to suppress the nodes for which the tem- 
perature is already known. 

We have the eguation 

■ ‘ ' ( 2 -«> 

The solution procedure, essentially iterative/ 
is as follows; 

The initial temperature field being specified, 
the value of conductivity matrix elements is found using 
the specified conductivity temperature relations . U sing 
these values of conductivity , the equation (2.46) is 
solved to obtain the temperature field at the end of the 
time interval 6t. The values accruing from the new 
temperature field are then the basis for calculation of 
conductivity values to be used in the computation of the 
temperature field during the succeeding iteration. The 
process is repeated for desired number of iterations or 
for tho desired time period. 

The matrix equation (2.46) can be rewritten to 
enable suppression of nodal equations with constant or 
fixed temperatures. Denoting the conductivity matrix 
by [a] and the heat capacity matrix by [b] we can 


write. 
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[a+b] = [b-a] ( 2 . 48 ) 

where?, solution for the vector T is sought and vector 
Tq is the temperature field, just computed or at the 


beginning of 

the iteration. 

If ■ T represents the fixed 

temperatures , 

('2 . 43 ) can be 

rewritten as. 


'^ll ®11 

^12 ®12 

1 T ) 

J 1 

) t 


^21 ”^2 1 

^22 ®22 

!f ) 


11 

®11 

^12 

®12 

'2 1 

+ B 2 ^ 

^2 2 

®2 2 


■ 

^ ) 
T \ 


f h) 



2 B ^2 


^®22 


j ^ 


which leads to, 


^‘^11 ® 11 ^ 1 '^1 


[a,, + B, J >, T, ■ 


11 "11-’ ( o) 

+ ? 

= - sCAjj] - [Aj^ + 

..... ( 2 . 49 ) 

In the solution of the above eguation it would 
be of considerably economy if use were made of the 
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symmetry oF matrices [a] and [b]. Also, these are 
banded matrices, enabling one to store elements within 
one half the band width alone entailing thereby savings 
in space and computational times. 

Use was made of the NAG (Numerical Algorithms 
Group) FORTRAN library [?6] routine F04ACF for the solu- 
tion of the system of equations (2.49). The triangular 
ring element mesh is shown in Fig. 2.7. 
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3.1 RESULTS AND DISCUSSIONS 

The diffusion model and its finite element 
formulation that have been developed in the previous 
chapter are used to find solutions for global and micro- 
temperature Fields in the half space X > 0, subject to 
a temperature-pulse boundary condition. System sensiti- 
vities to non-linearities arising from temperature depen- 
dent thermal properties have been accounted for. Two 
sets of materials/ Graphite/Epoxy and carbon/carbon com- 
posites have been considered as specific examples. A 
PEM computer package is developed for this purpose and 
executed on the DEC 1090 svstem. This is attached as 
Appendix 1. The results have been presented in the 
graphical form and are discussed below. The temperature 
dependence of thermal properties have been shown in 
Appendices 2 and 3. 

The average values of the nodal temperatures at 
each cross-section of the fibre and the matrix have been 
plotted egalnet the axial length of the Graph Ite/Epoxy 
composite. In Pig. 3.1. This has been done for a 
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particulrir instant of time just before the withdrawal of 
the constant temperature source at the boundarv. 
represents the fiber temperature and the matrix tem- 

perature. It is seen that the fiber is at a higher tem- 
perature than the matrix through a major length of the 
composite. It tends to the matrix or composite tempera- 
ture at distances away from the source. It can also be 
observed that the non-linearity of temperature properties 
provide lower temperature values In comparison to the 
linear case. This shows that the decrease in thermal 
conductivity with increase in temperature for Graphite 
plays a dominant role in the resulting temperature dis- 
tribution, 

‘ The average temperatures as explained above have 
also been plotted (Fig. 3.2) for the situation just after 
the withdrawal of the constant temperature source in order 
to study the decay of the temperature fields in the two 
materials of the composite. Since the thermal conduc- 
tivity of the matrix material is less than that of the 
fiber material, the temperature drop in the matrix mate- 
rial is much slower and consequently it is at higher 
temperature near the boundary. This also accounts for 
greater temperatures of the fiber at intermediate regions. 
It can be noticed that non-linearities provide for faster 
growth and decay of the temperature fields. 



45 


t 


bohaviour Of average temperatures and 

T " at an intormedlate oross-section with the passage of 

time is shown in Pin 

evident, temperature builds 

up and dooavs faster In the fiber than In the matrix 

accounting for the higher thermal conductivity of the 

fiber material. The figure also shows the difference 

between the linear and non-linear behaviour of the com- 

pos i to. 

The ability of the analysis to predict the tem- 
poraturo microstructure is illustrated in Fig. 3.4. The 
figure showvs the temperature fields across an intermediate 
cross-section a short while before and after the with- 
drawal of the constant temperature source. 

A short while before removal of the source, the 
temperature across the fiber is at a higher value than 
that cicross the matrix and is practically uniform. In 
the matrix core, it gradually reduces from a value of the 
fiber temperature at the interface to a minimum at the 
unit cell boundary. This is due to the lower thermal 
conductivity of the matrix material. 

A little while after removal of the source, 
however, the trend reverses. The fiber, having greater 
thermal conductivity, loses its heat faster and Is conse- 
quently at a lower temperature than the matrix. 
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The material properties used in deriving the 
results discussed above represent a Graph it e/Fpoxy fiber 
reinforced composite. In this case, the thermal proper- 
ties depend rathor weakly on temperature. The results 
however, imply that temperature dependence of thermal 
properties can significantly influence the results. In 
an effort to demonstrate this point, a similar initial 
Value problem was posed and solved for a carbon/carbon 
composite. The material properties used in the evalua- 
tion arc provided in appendix III. The following dis- 
cussions reflect the observations made. 

Figure 3.5 represents the temperature distribu- 
tion along the axial length of the composite before with- 
drawal of the source. The results indicate that non- 
linear iti os have significantly altered the distributions. 
It is worthwhile."’ to note that temperatures calculated by 
taking into consideration the non-linear material beha- 
viour, which is representative of actual working condi- 
tions to a greater extent, are much lower than those 
obtained by the assumption of linear behaviour. A possi- 
ble consequence of this knowledge would perhaps lead us 
to employ lesser insulating thicknesses than we possibly 
would have, entailing cost savings. 
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On withdrawal of the constant temperature source 
the tonipornture distributions are greatly altered as can 
be seen from the figure 3,6. The figure also provides 
ample evidence of the effect of non-linearities on the 
temperature distribution. It can be observed that 'hon- 
linear'^ temperatures are much lower than the " linear" 
temperatures for most of the length of the composite. 

Figure 3.7 shows the temperature time history at 
an intermedlnte cross-section of the composite. As long 
as tine .^'.ource exists, the temperatures build up. On 
removal of the source however, there is a decay in the 
temperature fields. The figure also iilusbtatos signifi- 
cant deviations caused due to non-linearities in the 
thermal bc?haviour of the composite. 

The temperature distribution across an inter- 
mediate scct.ion, just before, and after withdrawal of the 
constant temTierature source are shown in figure 3.8. As 
pxprertod, It is observed that the temperature builds up 
and df.'cavs faster across the cross— sect io-n of the fiber 
due to the higher thermal conductivity of the fiber mate- 
rial as compared to that of the matrix material. 
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3.2 CONCLUSIONS 

The beet conduction in uni-directional fiber 
reinforced composites subjected to a temperature pulse 
has been discussed in the present work. The effect of 
non-linearities on the temperature distributions have 
been amply demonstrated. The following conclusions may 
be drawn: 

(i) Non-1 Lnearities arising from temperature- dependent 
thermal properties significantly influence the 
temperature distributions. 

(ii) The non-linearities in case of both Granhite/ 

Epoxy and Carbon/Carbon, lead to temperature 
fields having magnitudes lower than those in 
the linear case.-. 

(iii) The non-linear analysis of the problem leads to 
an important result that less composite material 
is actually needed for insulation purposes than 
what is predicted by the linear solutions - 

(iv) Since the non-linear behaviour of the composite 
structures is a true representation of the 
materials used in practice, it is recommended 
that non-linearities must be accounted for in 
the heat diffusion analysis. 
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APPENDIX 2 


Material properties for Graphite/ Epoxy composite 
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APPENDIX 3 
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Material properties of carton/carbon composite 
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